{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topological feature extraction from graphs\n",
    "\n",
    "``giotto-tda`` can extract topological features from undirected or directed graphs represented as adjacency matrices, via the following transformers:\n",
    "\n",
    "- [VietorisRipsPersistence](https://giotto-ai.github.io/gtda-docs/latest/modules/generated/homology/gtda.homology.VietorisRipsPersistence.html) and [SparseRipsPersistence](https://giotto-ai.github.io/gtda-docs/latest/modules/generated/homology/gtda.homology.VietorisRipsPersistence.html) initialized with ``metric=\"precomputed\"``, for undirected graphs;\n",
    "- [FlagserPersistence](https://giotto-ai.github.io/gtda-docs/latest/modules/generated/homology/gtda.homology.VietorisRipsPersistence.html) initialized with ``directed=True``, for directed graphs, and with ``directed=False`` for undirected ones.\n",
    "\n",
    "In this notebook, we build some basic intuition on how these methods are able to capture structures and patterns from such graphs. We will focus on the multi-scale nature of the information contained in the final outputs (\"persistence diagrams\"), as well as on the differences between the undirected and directed cases. Although adjacency matrices of sparsely connected and even unweighted graphs can be passed directly to these transformers, they are interpreted as *weighted* adjacency matrices according to some non-standard conventions. We will highlight these conventions below.\n",
    "\n",
    "The mathematical technologies used under the hood are various flavours of \"persistent homology\" (as is also the case for [EuclideanCechPersistence](https://giotto-ai.github.io/gtda-docs/latest/modules/generated/homology/gtda.homology.EuclideanCechPersistence.html) and [CubicalPersistence](https://giotto-ai.github.io/gtda-docs/latest/modules/generated/homology/gtda.homology.CubicalPersistence.html)). If you are interested in the details, you can start from the [theory glossary](https://giotto-ai.github.io/gtda-docs/latest/theory/glossary.html) and references therein.\n",
    "\n",
    "If you are looking at a static version of this notebook and would like to run its contents, head over to [GitHub](https://github.com/giotto-ai/giotto-tda/blob/master/examples/persistent_homology_graphs.ipynb) and download the source.\n",
    "\n",
    "\n",
    "## See also\n",
    "\n",
    "- [Topological feature extraction using VietorisRipsPersistence and PersistenceEntropy](https://giotto-ai.github.io/gtda-docs/latest/notebooks/vietoris_rips_quickstart.html) which treats the \"special case\" of point clouds (see below).\n",
    "- [Plotting in giotto-tda](https://giotto-ai.github.io/gtda-docs/latest/notebooks/plotting_api.html), particularly Section 1.2 (as above, treats the case of point clouds).\n",
    "- [Classifying 3D shapes](https://giotto-ai.github.io/gtda-docs/latest/notebooks/classifying_shapes.html) (a more advanced example).\n",
    "- [Computing persistent homology of directed flag complexes](https://arxiv.org/abs/1906.10458) by Daniel Luetgehetmann, Dejan Govc, Jason Smith, and Ran Levi. \n",
    "\n",
    "**License: AGPLv3**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy.random import default_rng\n",
    "rng = default_rng(42)  # Create a random number generator\n",
    "\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "from scipy.sparse import coo_matrix\n",
    "\n",
    "from gtda.graphs import GraphGeodesicDistance\n",
    "from gtda.homology import VietorisRipsPersistence, SparseRipsPersistence, FlagserPersistence\n",
    "\n",
    "from igraph import Graph\n",
    "\n",
    "from IPython.display import SVG, display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Undirected graphs – ``VietorisRipsPersistence`` and ``SparseRipsPersistence``\n",
    "\n",
    "### General API\n",
    "\n",
    "If you have a collection ``X`` of adjacency matrices of graphs, you can instantiate transformers of class ``VietorisRipsPersistence`` or ``SparseRipsPersistence`` by setting the parameter ``metric`` as ``\"precomputed\"``, and then call ``fit_transform`` on ``X` to obtain the corresponding collection of *persistence diagrams* (see **Understanding the computation** below for an explanation).\n",
    "\n",
    "In the case of ``VietorisRipsPersistence``, ``X`` can be a list of sparse or dense matrices, and a basic example of topological feature extraction would look like this:\n",
    "```\n",
    "# Instantiate topological transformer\n",
    "VR = VietorisRipsPersistence(metric=\"precomputed\")\n",
    "\n",
    "# Compute persistence diagrams corresponding to each graph in X\n",
    "diagrams = VR.fit_transform(X)\n",
    "```\n",
    "\n",
    "Each entry in the result can be plotted as follows (where we plot the 0th entry, i.e. `diagrams[0]`):\n",
    "```\n",
    "VR.plot(diagrams, sample=0)\n",
    "```\n",
    "\n",
    "*Note*: ``SparseRipsPersistence`` implements an approximate scheme for computing the same topological quantities as ``VietorisRipsPersistence``. This can be useful for speeding up the computation on large inputs, but we will not demonstrate its use in this notebook.\n",
    "\n",
    "### Fully connected and weighted\n",
    "\n",
    "We now turn to the case of fully connected and weighted (FCW) undirected graphs. In this case, the input should be a list of 2D arrays or a single 3D array. Here is a simple example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a single weighted adjacency matrix of a FCW graph\n",
    "n_vertices = 10\n",
    "x = rng.random((n_vertices, n_vertices))\n",
    "# Fill the diagonal with zeros (not always necessary, see below)\n",
    "np.fill_diagonal(x, 0)\n",
    "\n",
    "# Create a trivial collection of weighted adjacency matrices, containing x only\n",
    "X = [x]\n",
    "\n",
    "# Instantiate topological transformer\n",
    "VR = VietorisRipsPersistence(metric=\"precomputed\")\n",
    "\n",
    "# Compute persistence diagrams corresponding to each entry (only one here) in X\n",
    "diagrams = VR.fit_transform(X)\n",
    "\n",
    "print(f\"diagrams.shape: {diagrams.shape} ({diagrams.shape[1]} topological features)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-fully connected weighted graphs\n",
    "\n",
    "In ``giotto-tda``, a non-fully connected weighted graph can be represented by an adjacency matrix in one of two possible forms:\n",
    "- a dense square array with ``np.inf`` in position $ij$ if the edge between vertex $i$ and vertex $j$ is absent.\n",
    "- a sparse matrix in which the non-stored edge weights represent absent edges.\n",
    "\n",
    "**Important notes**\n",
    "- A ``0`` in a dense array, or an explicitly stored ``0`` in a sparse matrix, does *not* denote an absent edge. It denotes an edge with weight 0, which, in a sense, means the complete opposite! See the section ***Understanding the computation*** below.\n",
    "- Dense Boolean arrays are first converted to numerical ones and then interpreted as adjacency matrices of FCW graphs. ``False`` values therefore should not be used to represent absent edges.\n",
    "\n",
    "### Understanding the computation \n",
    "\n",
    "To understand what these persistence diagrams are telling us about the input weighted graphs, we briefly explain the **clique complex (or flag complex) filtration** procedure underlying the computations in ``VietorisRipsPersistence`` when ``metric=\"precomputed\"``, via an example.\n",
    "\n",
    "Let us start with a special case of a weighted graph with adjacency matrix as follows:\n",
    "\n",
    "- the diagonal entries (\"vertex weights\") are all zero;\n",
    "- all off-diagonal entries (edge weights) are non-negative;\n",
    "- some edge weights are infinite (or very very large).\n",
    "\n",
    "We can lay such a graph on the plane to visualise it, drawing only the finite edges:\n",
    "\n",
    "![A weighted graph](images/weighted_graph.png)\n",
    "\n",
    "The procedure can be explained as follows: we let a parameter $\\varepsilon$ start at 0, and as we increase it all the way to infinity we keep considering the instantaneous subgraphs made of a) all the vertices in the original graph, and b) those edges whose weight is less than or equal to the current $\\varepsilon$. We also promote these subgraphs to more general structures called **(simplicial) complexes** that, alongside vertices and edges, also possess $k$**-simplices**, i.e. selected subsets of $k + 1$ vertices (a 2-simplex is an abstract \"triangle\", a 3-simplex an abstract \"tetrahedron\", etc). Our criterion is this: for each integer $k \\geq 2$, all $(k + 1)$-cliques in each instantaneous subgraph are declared to be the $k$-simplices of the subgraph's associated complex. By definition, the $0$-simplices are the vertices and the $1$-simplices are the available edges.\n",
    "\n",
    "As $\\varepsilon$ increases from 0 (included) to infinity, we record the following information:\n",
    "\n",
    "1. How many new **connected components** are created because of the appearance of vertices (in this example, all vertices \"appear\" in one go at $\\varepsilon = 0$, by definition!), or merge because of the appearance of new edges.\n",
    "2. How many new 1-dimensional \"holes\", 2-dimensional \"cavities\", or more generally $d$-dimensional **voids** are created in the instantaneous complex. A hole, cavity, or $d$-dimensional void is such only if there is no collection of \"triangles\", \"tetrahedra\", or $(d + 1)$-simplices which the void is the \"boundary\" of. *Note*: Although the edges of a triangle *alone* \"surround a hole\", these cannot occur in our particular construction because the \"filling\" triangle is also declared present in the complex when all its edges are.\n",
    "3. How many $d$-dimensional voids which were present at earlier values of $\\epsilon$ are \"filled\" by $(d + 1)$-simplices which just appear.\n",
    "\n",
    "This process of recording the full topological history of the graph across all edge weights is called (Vietoris-Rips) **persistent homology**.\n",
    "\n",
    "Let us start at $\\varepsilon = 0$: Some edges had zero weight in our graph, so they already appear!\n",
    "\n",
    "![$\\varepsilon = 0$](images/clique_complex_0_small.png)\n",
    "\n",
    "There are 9 connected components, and nothing much else.\n",
    "\n",
    "Up to and including $\\varepsilon = 2$, a few more edges are added which make some of the connected components merge together but do not create any hole, triangles, or higher cliques. Let us look at $\\varepsilon = 3$:\n",
    "\n",
    "![$\\varepsilon = 3$](images/clique_complex_3_small.png)\n",
    "\n",
    "The newly arrived edges reduce the number of connected components further, but more interestingly they create a 1D hole!\n",
    "\n",
    "As an example of a \"higher\"-simplex, at $\\varepsilon = 4$ we get our first triangle:\n",
    "\n",
    "![$\\varepsilon = 4$](images/clique_complex_4_small.png)\n",
    "\n",
    "At $\\varepsilon = 5$, our 1D hole is filled:\n",
    "\n",
    "![$\\varepsilon = 5$](images/clique_complex_5_small.png)\n",
    "\n",
    "At $\\varepsilon = 8$, two new 1D holes appear:\n",
    "\n",
    "![$\\varepsilon = 8$](images/clique_complex_8_small.png)\n",
    "\n",
    "Finally, at $\\varepsilon = 9$, some more connected components merge, but no new voids are either created or destroyed:\n",
    "\n",
    "![$\\varepsilon = 9$](images/clique_complex_9_small.png)\n",
    "\n",
    "We can stop as we have reached the maximum value of $\\varepsilon$, beyond which nothing will change: there is only one connected component left, but there are also two 1D holes which will never get filled.\n",
    "\n",
    "Fit-transforming via ``VietorisRipsPersistence(metric=\"precomputed\")`` on the original graph's adjacency matrix would return the following 3D array of **persistence diagrams**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diagrams = np.array([[[0., 1., 0],\n",
    "                      [0., 2., 0],\n",
    "                      [0., 2., 0],\n",
    "                      [0., 3., 0],\n",
    "                      [0., 4., 0],\n",
    "                      [0., 5., 0],\n",
    "                      [0., 6., 0],\n",
    "                      [0., 7., 0],\n",
    "                      [3., 5., 1],\n",
    "                      [8., np.inf, 1],\n",
    "                      [8., np.inf, 1]]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The notebook [Topological feature extraction using VietorisRipsPersistence and PersistenceEntropy](https://giotto-ai.github.io/gtda-docs/latest/notebooks/vietoris_rips_quickstart.html) explains how to interpret this output and how to make informative 2D scatter plots out of its entries. Here, we only have one entry corresponding to our graph:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gtda.plotting import plot_diagram\n",
    "\n",
    "plot_diagram(diagrams[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Small aside*: You would be correct to expect an additional row ``[0, np.inf, 0]`` representing one connected component which lives forever. By convention, since such a row would always be present under this construction and hence give no useful information, all transformers discussed in this notebook remove this feature from the output.\n",
    "\n",
    "#### Advanced discussion: Non-zero vertex weights and negative edge weights\n",
    "\n",
    "Although we introduced the simplifying assumptions that the diagonal entries of the input adjacency matrix is zero, and that all edge weights are non-negative, for the procedure to make sense we need a lot less. Namely:\n",
    "- The diagonal entry corresponding to a vertex is always interpreted as the value of the parameter $\\varepsilon$ at which that vertex \"appears\". Making all these entries equal to zero means, as in the example above, that all vertices appear simultaneously at $\\varepsilon = 0$. Generally however, different vertices can be declared to \"appear\" at different values, and even at negative ones.\n",
    "- The only constraint on each edge weight is that it should be no less than the \"vertex weight\" of either of its boundary vertices.\n",
    "\n",
    "As a simple example, subtracting a constant from *all* entries of an adjacency matrix has the effect of shifting all birth and death values by the same constant.\n",
    "\n",
    "### The \"special case\" of point clouds\n",
    "\n",
    "The use of ``VietorisRipsPersistence`` to compute multi-scale topological features of concrete point clouds in Euclidean space is covered briefly in Section 1.2 of [Plotting in giotto-tda](https://giotto-ai.github.io/gtda-docs/latest/notebooks/plotting_api.html), and more in-depth in [Classifying 3D shapes](https://giotto-ai.github.io/gtda-docs/latest/notebooks/classifying_shapes.html) and in [Can two-dimensional topological voids exist in two dimensions?](https://giotto-ai.github.io/gtda-docs/latest/notebooks/voids_on_the_plane.html)\n",
    "\n",
    "The Vietoris-Rips procedure for point clouds is often depicted as a process of growing balls of ever increasing radius $r$ around each point, and drawing edges between two points whenever their two respective $r$-balls touch for the first time. Just as in our clique complex construction above, cliques present at radius $r$ are declared to be higher-dimensional simplices in the instantaneous complex:\n",
    "\n",
    "![Vietoris–Rips filtration of a point cloud](images/vietoris_rips_point_cloud.gif)\n",
    "\n",
    "And just as in the case of weighted graphs, we record the appearance/disappearance of connected components and voids as we keep increasing $r$.\n",
    "\n",
    "The case of point clouds can actually be thought of as a special case of the case of FCW graphs. Namely, if:\n",
    "\n",
    "1. we regard each point in the cloud as an abstract vertex in a graph,\n",
    "2. we compute the square matrix of pairwise (Euclidean or other) distances between points in the cloud, and\n",
    "3. we run the procedure explained above with $\\varepsilon$ defined as $2r$,\n",
    "then we compute exactly the \"topological summary\" of the point cloud.\n",
    "\n",
    "So, in ``giotto-tda``, we can do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 1 point cloud with 20 points in 5 dimensions\n",
    "point_cloud = rng.random((20, 5))\n",
    "# Corresponding matrix of Euclidean pairwise distances\n",
    "pairwise_distances = squareform(pdist(point_cloud))\n",
    "\n",
    "# Default parameter for ``metric`` is \"euclidean\"\n",
    "X_vr_pc = VietorisRipsPersistence().fit_transform([point_cloud])\n",
    "\n",
    "X_vr_graph = VietorisRipsPersistence(metric=\"precomputed\").fit_transform([pairwise_distances])\n",
    "\n",
    "assert np.array_equal(X_vr_pc, X_vr_graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unweighted graphs and chaining with ``GraphGeodesicDistance``\n",
    "\n",
    "What if, as is the case in many applications, our graphs have sparse connections and are unweighted?\n",
    "\n",
    "In ``giotto-tda``, there are two possibilities:\n",
    "\n",
    "1. Encode the graphs as adjacency matrices of non-fully connected weighted graphs, where all weights corresponding to edges which are present are equal to ``1.`` (or any other positive constant). See section ***Non-fully connected weighted graphs*** above for the different encoding conventions for sparse and dense matrices.\n",
    "2. Preprocess the unweighted graph via [GraphGeodesicDistance](https://giotto-ai.github.io/gtda-docs/latest/modules/generated/graphs/processing/gtda.graphs.GraphGeodesicDistance.html) to obtain a FCW graph where edge $ij$ has as weight the length of the shortest path from vertex $i$ to vertex $j$ (and ``np.inf`` if no path exists between the two vertices in the original graph).\n",
    "\n",
    "### Example 1: Circle graph\n",
    "\n",
    "We now explore the difference between the two approaches in the simple example of a circle graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function -- directed circles will be needed later\n",
    "def make_circle_adjacency(n_vertices, directed=False):\n",
    "    weights = np.ones(n_vertices)\n",
    "    rows = np.arange(n_vertices)\n",
    "    columns = np.arange(1, n_vertices + 1) % n_vertices\n",
    "    return coo_matrix((weights, (rows, columns)))\n",
    "\n",
    "n_vertices = 10\n",
    "undirected_circle = make_circle_adjacency(n_vertices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can produce an SVG of the circle using ``python-igraph``, and display it.\n",
    "\n",
    "*Note*: If running from a live jupyter session, this will dump a file inside your notebook's directory. If ``pycairo`` is installed, you can draw the graph directly in the notebook by instead running\n",
    "```\n",
    "from igraph import plot\n",
    "plot(graph)\n",
    "```\n",
    "in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row, col = undirected_circle.nonzero()\n",
    "graph = Graph(n=n_vertices, edges=list(zip(row, col)), directed=False)\n",
    "fname = \"undirected_circle.svg\"\n",
    "graph.write_svg(fname)\n",
    "display(SVG(filename=fname))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Approach 1 means passing the graph as is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VietorisRipsPersistence(metric=\"precomputed\").fit_transform_plot([undirected_circle]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circular nature has been captured by the single point in homology dimension 1 ($H_1$) which is born at 1 and lives forever.\n",
    "\n",
    "Compare with what we observe when preprocessing first with ``GraphGeodesicDistance``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ggd = GraphGeodesicDistance(directed=False, unweighted=True).fit_transform([undirected_circle])\n",
    "VietorisRipsPersistence(metric=\"precomputed\").fit_transform_plot(X_ggd);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is still a \"long-lived\" topological feature in dimension 1, but this time its death value is finite. This is because, at some point, we have enough triangles to completely fill the 1D hole. Indeed, when the number of vertices/edges in the circle is large, the death value is around one third of this number. So, relative to the procedure without ``GraphGeodesicDistance``, the death value now gives additional information about the *size* of the circle graph!\n",
    "\n",
    "### Example 2: Two disconnected circles\n",
    "\n",
    "Suppose our graph contains two disconnected circles of different sizes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_vertices_small, n_vertices_large = n_vertices, 2 * n_vertices\n",
    "undirected_circle_small = make_circle_adjacency(n_vertices_small)\n",
    "undirected_circle_large = make_circle_adjacency(n_vertices_large)\n",
    "row_small, col_small = undirected_circle_small.nonzero()\n",
    "row_large, col_large = undirected_circle_large.nonzero()\n",
    "row = np.concatenate([row_small, row_large + n_vertices])\n",
    "col = np.concatenate([col_small, col_large + n_vertices])\n",
    "data = np.concatenate([undirected_circle_small.data, undirected_circle_large.data])\n",
    "two_undirected_circles = coo_matrix((data, (row, col)))\n",
    "\n",
    "graph = Graph(n=n_vertices_small + n_vertices_large, edges=list(zip(row, col)), directed=False)\n",
    "fname = \"two_undirected_circles.svg\"\n",
    "graph.write_svg(fname)\n",
    "display(SVG(filename=fname))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, the first procedure just says \"there are two 1D holes\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VietorisRipsPersistence(metric=\"precomputed\").fit_transform_plot([two_undirected_circles]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second procedure is again much more informative, yielding a persistence diagram with two points in homology dimension 1 with different finite deaths, each corresponding to one of the two circles:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ggd = GraphGeodesicDistance(directed=False, unweighted=True).fit_transform([two_undirected_circles])\n",
    "VietorisRipsPersistence(metric=\"precomputed\").fit_transform_plot(X_ggd);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Directed graphs – ``FlagserPersistence``\n",
    "\n",
    "Together with the companion package ``pyflagser`` ([source code](https://github.com/giotto-ai/pyflagser), [API reference](https://docs-pyflagser.giotto.ai/)), ``giotto-tda`` can extract topological features from *directed* graphs via the ``FlagserPersistence`` transformer.\n",
    "\n",
    "Unlike ``VietorisRipsPersistence`` and ``SparseRipsPersistence``, ``FlagserPersistence`` *only* works on graph data, so there is no ``metric`` parameter to be set. The conventions on input data are the same as in the undirected case, cf. section ***Non-fully connected weighted graphs*** above.\n",
    "\n",
    "The ideas and constructions underlying the algorithm in this case are very similar to the ones described above for the undirected case. Again, we threshold the graph and its directed edges according to an ever-increasing parameter and the edge weights. And again we look at \"cliques\" of vertices to define simplices and hence a \"complex\" for each value of the parameter. The main difference is that here simplices are **ordered** sets (tuples) of vertices, and that in each instantaneous complex the \"clique\" $(v_0, v_1, \\ldots, v_k)$ is a $k$-simplex if and only if, for each $i < j$, $(v_i, v_j)$ is a currently present directed edge.\n",
    "\n",
    "| (1, 2, 3) ***is not*** a 2-simplex in the complex | (1, 2, 3) ***is*** a 2-simplex in the complex |\n",
    "| :-- | :-- |\n",
    "| ![Directed flag complex with a hole](images/nontrivial_cycle_directed_flag_complex.svg) | ![Directed flag complex without a hole](images/simplex_directed_flag_complex.svg) |\n",
    "| (1, 2), (2, 3) and (3, 1) **form a 1D hole** | (1, 2), (2, 3) and (1, 3) form the boundary of (1, 2, 3) – **not a 1D hole** |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has interesting consequences: in the examples above, the left complex, in which the edges of the triangle \"loop around\" in the same direction, contains a 1D hole. On the other hand, the right one does not!\n",
    "\n",
    "### Example 1: Directed circle\n",
    "\n",
    "Let's try this on a \"directed\" version of the circle from earlier:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_vertices = 10\n",
    "\n",
    "directed_circle = make_circle_adjacency(n_vertices, directed=True)\n",
    "row, col = directed_circle.nonzero()\n",
    "\n",
    "graph = Graph(n=n_vertices, edges=list(zip(row, col)), directed=True)\n",
    "fname = \"directed_circle.svg\"\n",
    "graph.write_svg(fname)\n",
    "display(SVG(filename=fname))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Passing this directly to ``FlagserPersistence`` gives an unsurprising result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FlagserPersistence().fit_transform_plot([directed_circle]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can chain with an instance of ``GraphGeodesicDistance`` to get more information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ggd = GraphGeodesicDistance(directed=True, unweighted=True).fit_transform([directed_circle])\n",
    "FlagserPersistence().fit_transform_plot(X_ggd);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that this time the death time of the circular feature is circa one half of the number of vertices/edges. Compare this with the one-third factor we observed in the case of ``VietorisRipsPersistence``.\n",
    "\n",
    "### Example 2: Circle with alternating edge directions\n",
    "\n",
    "What happens when we make some of the edges flow the other way around the circle?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_flipped = np.concatenate([row[::2], col[1::2]])\n",
    "column_flipped = np.concatenate([col[::2], row[1::2]])\n",
    "\n",
    "graph = Graph(n=n_vertices, edges=list(zip(row_flipped, column_flipped)), directed=True)\n",
    "fname = \"directed_circle.svg\"\n",
    "graph.write_svg(fname)\n",
    "display(SVG(filename=fname))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the adjacency matrix\n",
    "weights = np.ones(n_vertices)\n",
    "directed_circle_flipped = coo_matrix((weights, (row_flipped, column_flipped)))\n",
    "\n",
    "# Run FlagserPersistence directly on the adjacency matrix\n",
    "FlagserPersistence().fit_transform_plot([directed_circle_flipped]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is identical to the persistence diagram for the directed circle (and for the undirected circle using ``VietorisRipsPersistence``). We cannot tell the difference between the two directed graphs when the adjacency matrices are fed directly to ``FlagserPersistence``. Let's try with ``GraphGeodesicDistance``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ggd = GraphGeodesicDistance(directed=True, unweighted=True).fit_transform([directed_circle_flipped])\n",
    "FlagserPersistence().fit_transform_plot(X_ggd);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in the case of the directed circle, the one-dimensional feature is born at 1. However, unlike that case, it persists all the way to infinity even after preprocessing with ``GraphGeodesicDistance``!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 3: Two oppositely-directed semicircles\n",
    "\n",
    "Our final example consists of a circle one half of which \"flows\" in one direction, and the other half in the other."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "row_two_semicircles = np.concatenate([row[:n_vertices // 2], col[n_vertices // 2:]])\n",
    "column_two_semicircles = np.concatenate([col[:n_vertices // 2], row[n_vertices // 2:]])\n",
    "\n",
    "graph = Graph(n=n_vertices, edges=list(zip(row_two_semicircles, column_two_semicircles)), directed=True)\n",
    "fname = \"two_directed_semicircles.svg\"\n",
    "graph.write_svg(fname)\n",
    "display(SVG(filename=fname))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the adjacency matrix\n",
    "weights = np.ones(n_vertices)\n",
    "two_semicircles = coo_matrix((weights, (row_two_semicircles, column_two_semicircles)))\n",
    "\n",
    "# Run FlagserPersistence directly on the adjacency matrix\n",
    "FlagserPersistence().fit_transform_plot([two_semicircles]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we passed the adjacency matrix directly and obtained the same persistence diagram as for the undirected circle. Let's try preprocessing with ``GraphGeodesicDistance``:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_ggd = GraphGeodesicDistance(directed=True, unweighted=True).fit_transform([two_semicircles])\n",
    "FlagserPersistence(directed=True).fit_transform_plot(X_ggd);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is similar to the persistence diagram for the coherently directed circle, but the death time for the topological feature in dimension 1 is slightly lower."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where to next?\n",
    "\n",
    "- Persistence diagrams are great for data exploration, but to feed their content to machine learning algorithms one must make sure the algorithm used is **independent of the relative ordering** of the birth-death pairs in each homology dimension. [gtda.diagrams](https://giotto-ai.github.io/gtda-docs/latest/modules/diagrams.html) contains a suite of vector representations, feature extraction methods and kernel methods that convert persistence diagrams into data structures ready for machine learning algorithms. Simple examples of their use are contained in [Topological feature extraction using VietorisRipsPersistence and PersistenceEntropy](https://giotto-ai.github.io/gtda-docs/latest/notebooks/vietoris_rips_quickstart.html), [Classifying 3D shapes](https://giotto-ai.github.io/gtda-docs/latest/notebooks/classifying_shapes.html) and [Lorenz attractor](https://giotto-ai.github.io/gtda-docs/latest/notebooks/lorenz_attractor.html).\n",
    "- In addition to ``GraphGeodesicDistance``, [gtda.graphs](https://giotto-ai.github.io/gtda-docs/latest/modules/graphs.html) also contains transformers for the creation of graphs from point cloud or time series data.\n",
    "- Despite the name, [gtda.point_clouds](https://giotto-ai.github.io/gtda-docs/latest/modules/point_clouds.html) contains transformers for the alteration of distance matrices (which are just adjacency matrices of weighted graphs) as a preprocessing step for persistent homology.\n",
    "- ``VietorisRipsPersistence`` builds on the [giotto-ph](https://github.com/giotto-ai/giotto-ph) and [ripser.py](https://ripser.scikit-tda.org/index.html) projects. The latter's website contains two tutorials on additional ways in which graphs can be constructed from [time series data](https://ripser.scikit-tda.org/notebooks/Lower%20Star%20Time%20Series.html) or [image data](https://ripser.scikit-tda.org/notebooks/Lower%20Star%20Image%20Filtrations.html), and fed to the clique complex filtration construction. With a few simple modifications, the code can be adapted to the API of ``VietorisRipsPersistence``."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
